In den Geowissenschaften ist die Fernerkundung die einzige Messtechnik, die eine vollständige meßtechnische Abdeckung auch globaler Skalen ermöglicht. Zur erfolgreichen Nutzung gehören die souverände Anwendung von existierenden und notwendigerweise auch die Anpassung dieser bzw.die Entwicklung eigener Methoden.

Einleitung

Ein wichtiger Anwendungsfall in der Geo- oder Umweltinformatik ist die Erfassung von Veränderungen mittels Satelliten-, Flugzeug- und Drohnenbildern (Change Detectin Analysis). Oft wird dies in Verbindung mit bio- und geophysikalischen oder auch vom Menschen verursachten Prozessen genutzt, um ein tieferes Verständnis und die Möglichkeit zur Entwicklung von Vorhersagemodellen zu erhalten. Dabei sind Bildanalysemethoden von hervorgehobener Bedeutung, um Informationen zu extrahieren, die es erlauben, die zugrunde liegenden Prozesse zu identifizieren. Ein ständig wachsender Bereich ist die Integration von immensen Datenmengen sog. Big Data in die Analysen.

Die nahezu exponetiell wachsende und bereits gegenwärtig kaum noch zu bewältigende Flut an verfügbaren Bild- und Fernerkundungsdaten muss einfach und automatisiert zugänglich sein um sowohl für den wissenschaftlichen Erkenntnisgewinn als auch für gesellschaftliche Zukunftsaufgaben nutzbar zu sein. Wir beginnen mit einer anwendungs orientierten Aufgabe in dieser Übung.

Informationen aus Bilddaten

Unverarbeitete Satellitenbilder sind nicht zwingend informativ. Unser Auge kann zwar ein Echtfarbenbild relativ schlüssig und intuitiv interpretieren, aber eine zuverlässige und reproduzierbare wissenschaftliche Interpretation erfordert andere Ansätze. Zudem können Bild analysemethoden zusätzliche,unsichtbare und spezifischere Informationen ableiten. Wir haben bereits einfache Indizes berechnet. Auch die Ableitung des NDVI oder der Oberflächenalbedo ist eine physikalisch begründete Umwandlung von Bildsignalen in eine Meßgröße.

Um nützliche oder zielführende Informationen, z. B. über die Bodenbedeckung in einem Gebiet, zu erhalten, müssen wir die Daten daher fragestellungszentriert analysieren. Der wohl bekannteste und gängigste Ansatz ist die überwachte Klassifizierung von Bilddaten in Kategorien die von Interesse sind.

Diese Übung führt Sie in die Klassifizierung von Satelliten- und Luftvermessungsdaten in R ein.

Wir werden uns mit den folgenden Themen beschäftigen:

Vorbereiten der Arbeitsumgebung und Laden der Daten
Digtalisierung von Trainingsbereichen
Unüberwachte/überwachte Klassifizierung

k-means (über RStoolbox)
Rekursive Partitionierung und Regressionsbäume (über caret)
Miximum Likelihood (über RStoolbox)

Kurze Einführung in Klassifkation

Überwachte Klassifizierung

Bei der überwachten Klassifizierung von Landbedeckungen wird aus einer begrenzten Menge Trainingsdaten zur Landbedeckung ein Modell abgeleitet, das die jeweilige Landbedeckung für den gesamten Datensatz vorhersagt. Die Landbedeckungstypen werden also a priori definiert, und das Modell versucht, diese Typen auf der Grundlage der Ähnlichkeit zwischen den Eigenschaften der Trainingsdaten und dem Rest des Datensatzes vorherzusagen.

Solche Klassifizierungen erfordern im Allgemeinen mindestens fünf Schritte:

Zusammenstellung eines umfassenden Eingabedatensatzes, der eine oder mehrere Rasterebenen enthält.
Auswahl von Trainingsgebieten, d.h. Teilmengen von Eingabedatensätzen, für die der Fernerkundungsexperte den Landbedeckungstyp kennt. Das Wissen über die Landbedeckung kann z.B. aus eigenen oder fremden in situ Beobachtungen, Managementinformationen oder anderen Fernerkundungsprodukten (z.B. hochauflösenden Luftbildern) gewonnen werden.
Training eines Modells unter Verwendung der Trainingsflächen. Zu Validierungszwecken werden die Trainingsflächen häufig in eine oder mehrere Test- und Trainingsstichproben unterteilt, um die Leistung des Modellalgorithmus zu bewerten.
Anwendung des trainierten Modells auf den gesamten Datensatz, d. h. Vorhersage der Bodenbedeckungsart auf der Grundlage der Ähnlichkeit der Daten an jedem Ort mit den Klasseneigenschaften des Trainingsdatensatzes.

Bitte beachten Sie, dass alle Arten der Klassifizierung eine gründliche Datenvor- und nachverarbeitung erfordern, die im Mittelpunkt der kommenden Kurseinheiten stehen wird.

Die folgende Abbildung zeigt die Schritte einer überwachten Klassifikation im Detail. Die optionalen Segmentierungsoperationen sind obligatorisch für objektorientierte Klassifizierungen, die sich nicht nur auf die Geometrie der Objekte, sondern auch auf die Werte jeder einzelnen Rasterzelle im Eingabedatensatz stützen. Um einzelne Objekte wie Häuser oder Baumkronen abzugrenzen, verwenden Fernerkundungsexperten Segmentierungsalgorithmen, die die Homogenität der Pixelwerte innerhalb ihrer räumlichen Nachbarschaft berücksichtigen.

Change Detection Waldveränderung Nord-West-Harz

In diesem Tutorium werden die Sentinel-2-Bilder aus der vorherigen Übung verwendet.

Schritt 1 - Einrichten des der Arbeitsumgebung

Sie können entweder die gespeicherten Daten aus der vorangegangenen Einheit verwenden oder einen neuen Raumausschnitt bestimmen, herunterladen und bearbeiten. Im Prinzip wird jedoch zuerst die Arbeitsumgebung geladen.

#  ---- 0 Projekt Setup ----
# Achtung Pakete müssen evtl. manuell installiert werden
library(envimaR)

#--- Schalter für den Download der sentinel daten
get_sen = FALSE

#--- schalter ob digitalisiert werden muss falls er auf FALSE gesetzt ist werden die
# (zuvor erstellten und gesciherten Daten ) im else Teil der Verzeigung eingelesen
digitize = FALSE

## setzen des aktuellen Projektverzeichnisses (erstellt mit envimaR) als root_folder
#root_folder = find_rstudio_root_file()
root_folder = "~/edu/geoinfo/"

# Einlesen des zuvor erstellten Setup-Skripts
source(file.path(root_folder, "src/functions/000_setup.R"))

nclasses=2

Bitte ergänzen Sie (bei auftreteneden Fehlermeldungen) etwaig fehlende oder defekte Pakete im obigen Setup-Skripts.

Auf der Grundlage der verfügbaren Sentinel Daten sollten nun als Erstes geeignete Datensätze für eine Oberflächenklassifikation identifiziert werden. Hierzu kann der vollständige Datensatz auch vom Kursdatenserver heruntergeladen werden (Bitte beachten Sie dass sie im VPN bzw. UniNetz angemeldet sein müssen). Entpacken Sie diese Daten in das Wurzelverzeichnis des Kursprojekts. Das heisst der dort bereits vorhandene Ordner data wird ersetzt/ergänzt.

Schritt 2 Übersicht verschaffen

Bei näherer Betrachtung der RGB Bilder (RGB432B) zeigt sich, das vier Datensätze aufgrund der Bildqualität und geringen Wolkenbedeckung geeignet zu sein scheinen. Es handelt sich um den 19.06. bzw. 24.7. 2019 und den 33.06 bzw. 30.7. 2020. Aufgrund des früheren Vegetationszeitpunktes wurden die Maibilder gewählt, da hier evtl. Aufwuchs auf den Rodungsflächen weniger stark sichtbar sein könnte.

Zunächst müssen diese Daten in einem “Rasterstapel” also einem Mehrkanalbild verfügbar sein. Bereits bei der ersten Beschäftigung mit dem sen2r Paket hatten wir weitere Produkte als sinnvoll erachtet und berechnet. Es sind folgende Indizes:

Normalized Difference Vegetation Index NDVI
Modified Soil Adjusted Vegetation Index MSAVI2
Simple Ratio 1600/820 Moisture Stress Index MSI
Modified Soil Adjusted Vegetation Index SAVI * Enhanced Vegetation Index EVI.

Diese Daten müssen für den gewünschten Zeitslot und Raumausschnitt besorgt und vorbereitet werden. Aus Gründen der Vereinfachung gibt es für den Kurs eine entsprechende sen2R-Projektdatei:

#--- Download der Daten
# gui = TRUE ruft die GUI zur Kontrolle auf
if (get_sen){
   out_paths_3 <- sen2r(
    gui = T,
    param_list = "~/edu/geoinfo/data/harz.json",
    tmpdir = envrmt$path_tmp,
  )
}
#--- Einlesen der Daten aus den Verzeichnissen
# RGB stack der beiden Jahre

pred_stack_2019 = raster::stack(list.files(file.path(envrmt$path_data_lev1,"RGB843B"),pattern = "20190619",full.names = TRUE))
pred_stack_2020 = raster::stack(list.files(file.path(envrmt$path_data_lev1,"RGB843B"),pattern = "20200623",full.names = TRUE))

# Stack-Loop über die Daten
for (pat in c("RGB432B","EVI","MSAVI2","NDVI","SAVI")){
  pred_stack_2019 = raster::stack(pred_stack_2019,raster::stack(list.files(file.path(envrmt$path_data_lev1,pat),pattern = "20190619",full.names = TRUE)))
  pred_stack_2020 = raster::stack(pred_stack_2020,raster::stack(list.files(file.path(envrmt$path_data_lev1,pat),pattern = "20200623",full.names = TRUE)))
}

# Zuweisen von leserlichen Namen auf die Datenebenen
names(pred_stack_2019) = c("nir","red","green","red","green","blue","EVI","MSAVI2","NDVI","SAVI")
names(pred_stack_2020) = c("nir","red","green","red","green","blue","EVI","MSAVI2","NDVI","SAVI")
saveRDS(pred_stack_2019,paste0(envrmt$path_data,"pred_stack_2019.rds"))
saveRDS(pred_stack_2020,paste0(envrmt$path_data,"pred_stack_2020.rds"))
pred_stack_2019 = readRDS(paste0(envrmt$path_data,"pred_stack_2019.rds"))
pred_stack_2020 = readRDS(paste0(envrmt$path_data,"pred_stack_2020.rds"))
names(pred_stack_2019) = c("nir","red_1","green_1","red_2","green_2","blue","EVI","MSAVI2","NDVI","SAVI")
names(pred_stack_2020) = c("nir","red_1","green_1","red_2","green_2","blue","EVI","MSAVI2","NDVI","SAVI")

# visuelle Überprüfung der stacks
plot(pred_stack_2019)

plot(pred_stack_2020)

# Einlesen der Corine Daten
# Für den Download ist ein Konto notwendig https://land.copernicus.eu/pan-european/corine-land-cover
# Daher die Daten manuell herunterladen und in das Verzeichnis kopieren und entpacken
# Dann auskommentierten Tail ausführen
# corine_eu = raster(file.path(envrmt$path_data_lev0,"u2018_clc2018_v2020_20u1_raster100m/DATA/U2018_CLC2018_V2020_20u1.tif"))
# tmp = projectRaster(pred_stack_2019[[1]],crs = crs(corine_eu))
# corine_crop = raster::crop(corine_eu,tmp)
# corine_utm = projectRaster(corine_crop,crs = crs(pred_stack_2019))
# corine = resample(corine_utm,pred_stack_2019[[1]])
# raster::writeRaster(corine,file.path(envrmt$path_data_lev0,"/corine.tif"),overwrite=TRUE)

# ansonsten den Beipieldatensatz laden
corine = raster::raster(file.path(envrmt$path_data_lev0,"corine.tif"))

# Erstellen einer Wald-Maske
# Agro-forestry areas code=22, Broad-leaved forest code=23,
# Coniferous forest code=24, Mixed forest code=25
mask = reclassify(corine,c(-100,22,0,22,26,1,26,500,0))

mapview(mask)

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 1693870 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  1693870 '

Schritt 3 erster Überblick mit einer K-Means Klassifikation

Die wohl bekannteste unüberwachte Klassifizierungstechnik ist das K-means-Clustering, das auch als der “einfachster Algorithmus des maschinellen Lernens” bezeichnet wird.

In unserem Beispiel (auf 2 Klassen angewendet und mit unsuperClass Funktion aus dem RStoolbox Paket ausgeführt) sieht das wie folgt aus:

## k-means über RStoolbox
# Modell
prediction_kmeans_2019 = unsuperClass(pred_stack_2019, nClasses = 2,norm = TRUE, algorithm = "MacQueen")
# Klassifikation
mapview(prediction_kmeans_2019$map, col = c('darkgreen', 'burlywood', 'green'))

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 1693870 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  1693870 '

prediction_kmeans_2020 = unsuperClass(pred_stack_2020, nClasses = 2,norm = TRUE, algorithm = "MacQueen")
mapview(prediction_kmeans_2020$map, col = c('darkgreen', 'burlywood', 'green'))

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 1693870 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  1693870 '

Schritt 4 Trainingsdaten erstellen

Der nächste Schritt ist hinischtlich seiner Durchführung in R optional, bietet aber die Möglichkeit, schnell und effektiv einige Trainingsflächen zu digitalisieren, ohne die RStudio-Welt zu verlassen. Für größere Aufgaben ist es unerlässlich, auf den hohen Komfort der QGIS-Arbeitsumgebung zurückzugreifen. Für diese Übung verwenden wir mapedit, ein kleines, aber feines Paket, das die Digitalisierung am Bildschirm in Rstudio oder in einem Browser ermöglicht. In Kombination mit mapview ist es auch sehr komfortabel für schnelle Digitalisierungsaufgaben. Besonders hilfreich ist die Möglichkeit, beliebige Farbkomposita als Grundlage einzulesen.

Farbkomposita für bessere Trainingsergebnisse

mapview::viewRGB(pred_stack_2019, r = 1, g = 2, b = 3, maxpixels =  1693870) +
  mapview::viewRGB(pred_stack_2019, r = 5, g = 6, b = 7, maxpixels =  1693870)

## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same

## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same

Verwenden Sie die Ebenensteuerung, um die Ebenen umzuschalten. Bei Echtfarbkompositen werden die sichtbaren Spektralkanäle Rot (B04), Grün (B03) und Blau (B02) den entsprechenden roten, grünen bzw. blauen Farbkanälen zugeordnet, wodurch ein quasi natürliches “farbiges” Bild der Oberfläche entsteht, wie es ein Mensch sehen würde, der auf dem Satelliten sitzt. Falschfarbenbilder werden häufig mit den Spektralkanälen für das nahe Infrarot, Rot und Grün erzeugt. Sie eignen sich hervorragend für die Einschätzung der Vegetation, da Pflanzen nahes Infrarot und grünes Licht reflektieren, während sie rotes Licht absorbieren (Red Edge Effect). Ein dichterer Pflanzenbewuchs ist dunkler rot. Städte und offener Boden sind grau oder hellbraun, und Wasser erscheint blau oder schwarz.

Quick & Dirty Trainingsdaten digitalisieren

Wir nehmen an, dass wir zwei Arten von Landbedeckung klassifizieren wollen: Abholzungen und Anderes. Jede Klasse wird einzeln digitalisiert. Im Anschluss werden in den Trainingsgebieten die Merkmale des jeweiligen Rasterstacks extrahiert und auf Fehlwerte bereiningt. Falls dieser Teil bereits absolviert wurde kann die logische Variable (am Anfang des Scripts definiert) digitize auf FALSE gesetzt werden und es wird dann der ELSE Teil der Verzweigung durchlaufen - also nur noch die existierenden Daten eingelesen.

#---- Digitalisierung der Trainingsdaten ----

if (digitize) {
  # Für die überwachte Klassifikation benötigen wir Trainingsgebiete. Sie können Sie wie nachfolgend digitalisieren oder alternativ z.B. QGis verwenden
  
  #--- 2019
  # Kahlschlag
  # Es wird das Falschfarbenkomosit in originaler Auflösung genutzt (maxpixels =  1693870)
  # Bitte beachten Sie dass es (1) deutlich länger lädt und (2) Vegetation in Rot dargestellt wird.
  # Die Kahlschäge sind jetzt grün
  train_area <- mapview::viewRGB(pred_stack_2019, r = 1, g = 2, b = 3, maxpixels =  1693870) %>% mapedit::editMap()
  # Hinzufügen der Attribute class (text) und id (integer)
  clearcut <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "clearcut", id = 1)
  
  # other: hier gilt es möglichst verteilt übers Bild möglichst alle nicht zu Kahlschlag  gehörenden Flächen zu erfassen.
  train_area <- mapview::viewRGB(pred_stack_2019, r = 1, g = 2, b = 3) %>% mapedit::editMap()
  other <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "other", id = 2)
  
  # rbind  kopiert die beiden obigen Vektorobjekte in eine Datei
  train_areas_2019 <- rbind(clearcut, other)
  
  # Umprojizieren auf die Raster Datei
  train_areas_2019 = sf::st_transform(train_areas_2019,crs = sf::st_crs(pred_stack_2019))
  
  
  # Extraktion der Trainingsdaten für die digitalisierten Flächen
  tDF_2019 = exactextractr::exact_extract(pred_stack_2019, train_areas_2019,  force_df = TRUE,
                                          include_cell = TRUE,include_xy = TRUE,full_colnames = TRUE,include_cols = "class")
  #  auch hier wieder zusamenkopieren in eine Datei
  tDF_2019 = dplyr::bind_rows(tDF_2019)
  
  # Löschen von etwaigen Zeilen die NA (no data) Werte enthalten
  tDF_2019 = tDF_2019[complete.cases(tDF_2019) ,]
  tDF_2019 = tDF_2019[ ,rowSums(is.na(tDF_2019)) == 0]
  
  # check der extrahierten Daten
  summary(tDF_2019)
  mapview(train_areas_2019)+pred_stack_2019[[1]]
  
  # Abspeichern als R-internes Datenformat
  # ist im Repo hinterlegt und kann desahlb (zeile drunter) eingeladen werden
  saveRDS(tDF_2019, paste0(envrmt$path_data,"train_areas_2019.rds"))
  
  
  # # ---- Das gleiche muss für 2020 wiederholt werden zum digitalisieren und extrahieren bitte ent-kommentieren ----
  
  # # Kahlschlag
  # train_area <- mapview::viewRGB(pred_stack_2020, r = 1, g = 2, b = 3,maxpixels =  1693870) %>% mapedit::editMap()
  # clearcut <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "clearcut", id = 1)
  # train_area <- mapview::viewRGB(pred_stack_2020, r = 1, g = 2, b = 3) %>% mapedit::editMap()
  # other <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "other", id = 2)
  # train_areas_2020 <- rbind(clearcut, other)
  # train_areas_2020 = sf::st_transform(train_areas_2020,crs = sf::st_crs(pred_stack_2020))
  # tDF_2020 = exactextractr::exact_extract(pred_stack_2020, train_areas_2020,  force_df = TRUE,
  #                                         include_cell = TRUE,include_xy = TRUE,full_colnames = TRUE,include_cols = "class")
  # tDF_2020 = dplyr::bind_rows(tDF_2020)
  # tDF_2020 = tDF_2020[  rowSums(is.na(tDF_2020)) == 0,]
  # saveRDS(tDF_2020, paste0(envrmt$path_data,"train_areas_2020.rds"))
  
} else {
  train_areas_2019 = readRDS(paste0(envrmt$path_data,"train_areas_2019.rds"))
  train_areas_2020 = readRDS(paste0(envrmt$path_data,"train_areas_2020.rds"))
  
}

Schritt 5 - Klassifizierung und Prüfen der Modellgüte

Maschinelles Lernen wie z. B. Random Forest, Algorithmen können so trainiert werden, dass sie automatisiert Muster in Daten finden. Solange die Trainingsdaten geeignet und repräsentativ sind , können auch Vorhersagen für Räume getroffen werden, für die keine Daten vorhanden sind.

Wir wollen nun die räumlichen Merkmale Kahlschlag/kein Wald mit Techniken des maschinellen Lernens vorhersagen und Methoden der Zufallsvalidierung anwenden.

Wir wollen ein Modell erstellen, das zwischen den Pixeln eines Kahlschlages und allen anderen Pixeln unterscheiden kann. Dazu nutzen wir die HEruntergeladnen Daten im Vorhersage Rasterstack.

Random forest

Random Forests können sowohl für Regressions- als auch für Klassifizierungsaufgaben verwendet werden, wobei letztere besonders in der Umwelt-Fernerkundung relevant sind. Wie bei jeder maschinellen Lernmethode lernt auch das Random-Forest-Modell, Muster und Strukturen in den Daten selbst zu erkennen.

Abbildung: Vereinfachte Darstellung der Klassifizierung von Daten durch Random Forest während des Trainings. Venkata Jagannath [CC BY-SA 4.0] via wikipedia.org

Ein Random-Forest-Algorithmus lernt über die Daten, indem er viele Entscheidungsbäume erstellt - daher auch der Name “Forest.” Für Klassifizierungsaufgaben nimmt der Algorithmus eine passende Instanz aus eines Baumes aus dem Trainingsdatensatz und weist dem Pixel dies Klasse zu. Dies wird mit allen verfügbaren Entschidungsbäuumen so gemacht und letzlich wird nach dem Winner takes it all Ansatz die Instanz der Klasse zugeordnet, die das Ergebnis der meisten Bäume ist.

Da der Random-Forest-Algorithmus Trainingsdaten benötigt, handelt es sich um eine überwachte Lernmethode. Das bedeutet, dass wir als Nutzer dem Algorithmus sagen müssen, was er vorhersagen soll. In unserem Fall müssen die Trainingsdaten verschiedene Kategorien oder Klassifizierungen der Bodenbedeckung (z. B. Kahlschlag, Rest) enthalten und mit mit einem Tag so gekennzeichnet sein.

## ---- Überwachte  mit Klassifikation Random Forest ----
  train_areas_2019 = readRDS(paste0(envrmt$path_data,"train_areas_2019.rds"))
  train_areas_2020 = readRDS(paste0(envrmt$path_data,"train_areas_2020.rds"))
## Hier wird der Random Forest über das Utility Paket caret aufgerufen
# Setzen eines "seed" ermöglicht reproduzierbaren Zufall
set.seed(123)

# Zufälliges Ziehen von 25% der Daten (Training/Test)
idx_2019 = createDataPartition(train_areas_2019$class,list = FALSE,p = 0.25)
trainDat_2019 = train_areas_2019[idx_2019,]
testDat_2019 = train_areas_2019[-idx_2019,]
idx_2020 = createDataPartition(train_areas_2020$class,list = FALSE,p = 0.25)
trainDat_2020 = train_areas_2020[idx_2020,]
testDat_2020 = train_areas_2020[-idx_2020,]

#  Response-Variable (=Spalte "class") muss den Datentyp "factor" haben
trainDat_2019$class <- as.factor(trainDat_2019$class)
testDat_2019$class <- as.factor(testDat_2019$class)
trainDat_2020$class <- as.factor(trainDat_2020$class)
testDat_2020$class <- as.factor(testDat_2020$class)


# Einstellungen Modelltraining: cross-validation, 10 Wiederholungen
ctrlh = trainControl(method = "cv",
                     number = 10,
                     savePredictions = TRUE)

#--- random forest model training
cv_model_2019 = train(trainDat_2019[,2:11], # in den Spalten 2 bis 20 stehen die Trainingsdaten (Prediktoren genannt)
                      trainDat_2019[,1],         # in der Spalte 1 steht die zu klassizierende Variable (Response genannt)
                      method = "rf",             # Methode hier rf für random forest
                      metric = "Kappa",          # Qualitäts/Performanzmaß KAppa
                      trControl = ctrlh,         # obig erzeugte Trainingssteuerung soll eingelsen werden
                      importance = TRUE)         # Die Bedeung der Variablen wird mit abgespeichert

cv_model_2019

## Random Forest 
## 
## 5996 samples
##   10 predictor
##    2 classes: 'clearcut', 'other' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 5397, 5396, 5397, 5396, 5396, 5396, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    2    0.9879911  0.8578546
##    6    0.9884903  0.8658856
##   10    0.9879903  0.8600513
## 
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 6.

Berechnung einer Wahrheitsmatrix zur Überprüfung der Modellgüte

Die Testdaten werden nun für die unabhängige Qualitätsprüfung des Modells verwendet. Eine Wahrheits- oder Konfusionsmatrix zeigt an, wie präzise das Modell die korrekten Klassen vorhersagt. Die Hauptdiagonale der Matrix zeigt die Fälle an, in denen das Modell korrekt ist. In diesem Falle mit nur 2 Klassen gilt der Sonderfall der Beurteilung eines binären Klassifikators. Speziell für die hier verwendete Funktion finden sich ausführliche Erläuertungen in der caret Hilfe. Die wesentlichen Ausagen sind :

‘Positive’ Class = clearcut: wird mit der Sensitivität (true positive rate) erfasst, die die Wahrscheinlichkeit angibt, mit der ein positives Objekt korrekt als positiv klassifiziert wird.
’ Negative Class’ = other: wird mit der Spezifität (true negative rate) erfasst und gibt die Wahrscheinlichkeit an, mit der ein negatives Objekt korrekt als negativ klassifiziert wird
Positive and negative predictive values geben für clearcut bzw. other an. Sie sind um die reale Häufikeit korrigiert und ein Schätzmass für die Präzision bzw. Pereformance des Modells. Trotz der hohen Werte sehen wir das die Klasse clearcut deutlich abfällt. Insgesamt kann das Modell jedoch als gut bezeichnet werden.

# ---- Berechnung der Konfusionsmatrix  ----
cm_rf <- confusionMatrix(data = predict(cv_model_2019, newdata = testDat_2019), testDat_2019$class)
cm_rf

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction clearcut other
##   clearcut      699   103
##   other         120 17063
##                                           
##                Accuracy : 0.9876          
##                  95% CI : (0.9859, 0.9892)
##     No Information Rate : 0.9545          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.8559          
##                                           
##  Mcnemar's Test P-Value : 0.284           
##                                           
##             Sensitivity : 0.85348         
##             Specificity : 0.99400         
##          Pos Pred Value : 0.87157         
##          Neg Pred Value : 0.99302         
##              Prevalence : 0.04554         
##          Detection Rate : 0.03887         
##    Detection Prevalence : 0.04459         
##       Balanced Accuracy : 0.92374         
##                                           
##        'Positive' Class : clearcut        
##

# Klassifikation (auch Vorhersage genannt)
prediction_rf_2019  = raster::predict(pred_stack_2019 ,cv_model_2019 )

#--- 2020
# Einstellungen wie 2019
cv_model_2020 = train(trainDat_2020[,2:11],
                      trainDat_2020[,1],
                      method = "rf",
                      metric = "Kappa",
                      trControl = ctrlh,
                      importance = TRUE)

cv_model_2020

## Random Forest 
## 
## 7845 samples
##   10 predictor
##    2 classes: 'clearcut', 'other' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 7061, 7060, 7060, 7060, 7061, 7060, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    2    0.9847046  0.9218005
##    6    0.9841944  0.9192613
##   10    0.9838122  0.9174582
## 
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

prediction_rf_2020  = raster::predict(pred_stack_2020 ,cv_model_2020)

Maximum Likelihood Klassifikation

Bei der Maximum-Likelihood-Klassifizierung wird davon ausgegangen, dass die Verteilung der Daten für jede Klasse und in jedem Kanal normal verteilt sind. Unter dieser Vorraussetzung wird die Wahrscheinlichkeit berechnet, dass ein bestimmtes Pixel zu einer bestimmten Klasse gehört. Da auch die Wahrscheinlichkeiten als Schwellenwert angegeben werden können werden ohne diese Einschränkung alle Pixel ungeachtet wie unwahrscheinlich zugeordnet. Jedes Pixel wird der Klasse zugeordnet, die die höchste Wahrscheinlichkeit aufweist (d. h. die maximale Wahrscheinlichkeit).

# ---- Maximum Likelihood Classification ----
# superClass() Funktion aus dem Paket RSToolbox umwandeln der Tabelle in das
# geforderte SpatialdataPoint Objekt Identifikation der Koordinaten
sp::coordinates(trainDat_2019) = ~X+Y
sp::coordinates(trainDat_2020) = ~X+Y
# Zuweisen der korrekten Projektion (hier aus dem zugehörigen Raster-Stack)
crs(trainDat_2019) = crs(pred_stack_2019)
crs(trainDat_2020) = crs(pred_stack_2020)

# superClass method "mlc" berechnet die Statisik und klassifiziert in einem aufruf
prediction_mlc_2019       <- superClass(pred_stack_2019, trainData = trainDat_2019[,1:11], responseCol = "class",
                                        model = "mlc", tuneLength = 1, trainPartition = 0.5,verbose = FALSE)

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

prediction_mlc_2020       <- superClass(pred_stack_2020, trainData = trainDat_2020[,1:11], responseCol = "class",
                                        model = "mlc", tuneLength = 1, trainPartition = 0.5,verbose = FALSE)

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

# Ergebnisse der rf Klassifikation
#prediction_mlc_2019$model
#prediction_mlc_2020$model


## ---- Visualisierung mit mapview ----
mapview::viewRGB(mask*pred_stack_2020, r = 4, g =5, b = 6,maxpixels =  1693870)+
  mapview(mask*prediction_rf_2019 , alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = TRUE) +
  mapview(mask*prediction_rf_2020, alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
  mapview(mask*prediction_mlc_2019$map,alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
  mapview(mask*prediction_mlc_2020$map,alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE)

## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same

Ressourcen

Betrachten sie sie nachfolgenden Ressourcen als Beispiele dafür, wie aus solchen Quellen (die alle mehr oder weniger technisch ähnlich sind) ein bestimmtes Set von Werkzeugen zur Bearbeitung einer Fragestellung “entsteht”. Dann nach viel Recherche und kritischer Prüfung wird ein aktueller Stand der Technik innerhalb der wissenschaftlcihen Gemeinschaft erreicht der sich zudem natürlich stetig weiterentwickelt.

Schauen sie einmal in die nachstehende Auswahl von Blogs und Hilfestellungen:

Digitalisierung von Trainingsdaten by the Forest Inventory and Remote Sensing department at the University of Goettingen (Germany)
Digitalisierung Turorial in the QGIS 3.16 documentation
Robert J. Hijmans rspatial - supervised classification
Ivan Lizarazo RPubs Tutorial
Sydney Goldstein blog
João Gonçalves supervised classification
Valentin Stefan pixel-based supervised classification

Die beiden ersten Beiträge sind rein technische Tutorials. Die weiteren Beiträge mischen technische Anleitungen mit konzeptionellen oder konkreten fachlichen Fragestellungen. Dennoch sind sie hier nicht als fachwissenschaftliche Referenz gedacht. Sie dienen vielmehr dazu zu zeigen wie technisches und konzeptionelles Verständnis schrittweise erarbeiet werden kann und unterstützen, durch “nachkochen” und anwenden, eigenständiger die Lösung von Fragestellungen anzugehen.

Ich möchte ausdrücklich Valentin Stefan den Autor des Blogbeitrags pixel-based supervised classification zitieren:

“[…] Betrachten Sie diesen Inhalt als einen Blogbeitrag und nichts weiter. Er erhebt nicht den Anspruch, eine erschöpfende Übung oder ein Ersatz für Ihr kritisches Denken zu sein. […]”